#include "Matriz.h"
#include "VectorMath.h"

//constructor
template<class T>
Matriz<T> :: Matriz(unsigned f, unsigned c){
    datos = new VectorMath<T> [f];
    for(unsigned i = 0; i < f; i++)
        datos[i].resize( c );
    filas = f;
    columnas = c;
}

//constructor de copia
template<class T>
Matriz<T> :: Matriz(const Matriz<T> &m){
    filas = m.filas;
    columnas = m.columnas;
    datos = new VectorMath<T> [filas];
    for(unsigned i = 0; i < filas; i++)
        datos[i].resize( columnas );

    for(unsigned i = 0; i < filas; i++)
        for(unsigned j = 0; j < columnas; j++)
            datos[i][j] = m.datos[i][j];
}

//destructor
template<class T>
Matriz<T> :: ~Matriz(){
    delete [] datos;
}

//getter filas
template<class T>
unsigned Matriz<T> :: getFilas(){
    return filas;
}

//getter columnas
template<class T>
unsigned Matriz<T> :: getColumnas(){
    return columnas;
}

//operator []
template<class T>
VectorMath<T> & Matriz<T> :: operator [](unsigned i){
    return datos[i];
}

//operator =
template<class T>
Matriz<T> & Matriz<T> :: operator = (const Matriz<T> &m){
    if(filas != m.filas || columnas != m.columnas){
        delete[] datos;

        datos = new VectorMath<T> [m.filas];
        for(unsigned i = 0; i < m.filas; i++)
            datos[i].resize( m.columnas );
        filas = m.filas;
        columnas = m.columnas;
    }

    for(unsigned i = 0; i < filas; i++)
        for(unsigned j = 0; j < columnas; j++)
            datos[i][j] = m.datos[i][j];
    return *this;
}

//operator +
template<class T>
Matriz<T> Matriz<T> :: operator + (const Matriz<T> &m) const{
    Matriz<T> aux(*this);
    if(m.filas == filas && m.columnas == columnas){
        for(unsigned i = 0; i < filas; i++)
            for(unsigned j = 0; j < columnas; j++)
                aux.datos[i][j] += m.datos[i][j];
    }
    return aux;
}

//operator -
template<class T>
Matriz<T> Matriz<T> :: operator - (const Matriz<T> &m) const{
    Matriz<T> aux(*this);
    if(m.filas == filas && m.columnas == columnas){
        for(unsigned i = 0; i < filas; i++)
            for(unsigned j = 0; j < columnas; j++)
                aux.datos[i][j] -= m.datos[i][j];
    }
    return aux;
}

//operator *
template<class T>
Matriz<T> Matriz<T> :: operator * (const Matriz<T> &m) const{
    T acum = 0;
    Matriz<T> aux(filas, m.columnas);
    
    if(columnas == m.filas){
            
        for(unsigned f = 0; f < filas; f++){
            for(unsigned c = 0; c < m.columnas; c++){
                acum = 0;
                for(unsigned a = 0; a < columnas; a++)
                    acum += datos[f][a] * m.datos[a][c];
                aux.datos[f][c] = acum;
            }
        }
    }
    return aux;
}

//operator * VectorMath
template<class T>
VectorMath<T> Matriz<T> :: operator *  (const VectorMath<T> &vec) const{
   if(columnas != vec.size()) exit(1);
           
    VectorMath<T> aux(filas);
    T suma;
           
    for(unsigned i = 0; i < filas; i++){
          suma = 0;
        for(unsigned j = 0; j < vec.size(); j++)
            suma += datos[i][j] * vec[j];
        aux[i] = suma;
    }
    
    return aux;
}

//operator ==
template<class T>
bool Matriz<T> :: operator == (const Matriz<T> &m) const{
    if(filas != m.filas || columnas != m.columnas) return false;
    for(unsigned i = 0; i < filas; i++)
        for(unsigned j = 0; j < columnas; j++)
            if(datos[i][j] != m.datos[i][j])
                return false;
    return true;
}     

//adjunta
template<class T>
Matriz<T> Matriz<T> :: adjunta(){
    if(filas != columnas){
        cout << "Solo las matrices cuadradas tienen adjunta.\n";
        exit(0);
    }

    int orden = filas;
    Matriz<T> adj(orden, orden);

    switch(orden){
        case 1: ; break; //la adjunta de una matriz de orden 1 es la misma matriz
        case 2: adj.datos[0][0] = datos[1][1];
                adj.datos[0][1] = - datos[1][0];
                adj.datos[1][0] = - datos[0][1];
                adj.datos[1][1] = datos[0][0];
                break;
        case 3: adj.datos[0][0] = datos[1][1] * datos[2][2] - datos[1][2] * datos[2][1];
                adj.datos[1][0] = datos[1][2] * datos[2][0] - datos[1][0] * datos[2][2];
                adj.datos[2][0] = datos[1][0] * datos[2][1] - datos[1][1] * datos[2][0];
                adj.datos[0][1] = datos[2][1] * datos[0][2] - datos[2][2] * datos[0][1];
                adj.datos[1][1] = datos[2][2] * datos[0][0] - datos[2][0] * datos[0][2];
                adj.datos[2][1] = datos[2][0] * datos[0][1] - datos[2][1] * datos[0][0];
                adj.datos[0][2] = datos[0][1] * datos[1][2] - datos[0][2] * datos[1][1];
                adj.datos[1][2] = datos[0][2] * datos[1][0] - datos[0][0] * datos[1][2];
                adj.datos[2][2] = datos[0][0] * datos[1][1] - datos[0][1] * datos[1][0];
                break;
        default: cout << "Calcular adjunta de orden > 3 es muy complejo";
    }
    return adj;
}

//determinante
template<class T>
T Matriz<T> :: determinante(){
    if(filas != columnas){
        cout << "Solo las matrices cuadradas tienen determinante.\n";
        return 0;
    }

    int orden = filas;
    T det = 0;
    
    switch(orden){
        case 1: det = datos[0][0]; break;
        case 2: det = datos[0][0] * datos[1][1] - datos[1][0] * datos[0][1];
                break;
        case 3: det += datos[0][0] * datos[1][1] * datos[2][2];
                det += datos[0][1] * datos[1][2] * datos[2][0];
                det += datos[0][2] * datos[1][0] * datos[2][1];
                det -= datos[0][0] * datos[1][2] * datos[2][1];
                det -= datos[0][1] * datos[1][0] * datos[2][2];
                det -= datos[0][2] * datos[1][1] * datos[2][0];
                break;
        default: cout << "Calcular determinante de orden > 3 es muy complejo";
    }
    return det;
}

//inversa
template<class T>
Matriz<T> Matriz<T> :: inversa(){
        
    if(filas == columnas){
        Matriz<T> matriz(  filas , columnas + 1 ),
                  inv( filas , columnas );
        
        for( unsigned ciclo = 0 ; ciclo < filas ; ciclo++ ){
                  
            for( unsigned i = 0 ; i < filas ; i++ )
                for(unsigned j = 0 ; j < columnas + 1 ; j++)
                    matriz[i][j] = (*this)[i][j];
                    
            for( unsigned i = 0 ; i < filas ; i++ )
                matriz[ i ][ columnas ] = 0;
                
            matriz[ ciclo ][ columnas ] = 1;
            
            matriz = matriz.escalonada();
            
            for( unsigned i = 0 ; i < filas ; i++ )
                inv[ i ][ ciclo ] = matriz[ i ][ columnas ];
            
        }
        
        return inv;
    
    }else
        cout << "filas != columnas";
        
    return (*this);
    
}

//escalonada Gauss-Jordan devuelve la identidad y las incognitas
template<class T>
Matriz<T> Matriz<T> :: escalonada(){
    Matriz<T> esc(*this);

    for( unsigned pivot = 0 ; pivot < filas ; pivot++ ){
        //Cuadrante 1 y 2
        if(pivot > 0){
            for( unsigned i = 0 ; i < pivot ; i++ ){
                for( unsigned j = 0 ; j < filas + 1 ; j++ ){
                    if( j != pivot )
                        esc[i][j] = esc[i][j] * esc[pivot][pivot] - 
                                                esc[i][pivot] * esc[pivot][j];
                }
            }
        }
        
        //Cuadrante 4
        for( unsigned i = pivot + 1 ; i < filas ; i++ ){
            for( unsigned j = pivot + 1 ; j < filas + 1 ; j++ ){
                    esc[i][j] = esc[i][j] * esc[pivot][pivot] - 
                                            esc[i][pivot] * esc[pivot][j];
            }
        }
        
        for( unsigned i = 0 ; i < filas ; i++ )
            if( i != pivot )
                esc[i][pivot] = 0;
    }

    for( unsigned pivot = 0 ; pivot < filas ; pivot++ ){
        esc[pivot][filas] = esc[pivot][filas] / esc[pivot][pivot];
        esc[pivot][pivot] = 1;
    }

    return esc;
}

//toString
template<class T>
string Matriz<T> :: toString(){
    string matToString;
    
    ostringstream stream;
    stream << filas << "\n";
    stream << columnas << "\n";
    
    for( unsigned i = 0 ; i < filas ; i++ )
        for( unsigned j = 0 ; j < filas ; j++ )
            stream << datos[i][j] << "\n";
    
    matToString = stream.str();

    return matToString;
}

//test
template<class T>
bool Matriz<T> :: test(){
    Matriz<float> mat(2,2);
    
    mat[0][0] = 1;
    mat[0][1] = 2;
    mat[1][0] = 3;
    mat[1][1] = 4;
    
    cout << "Matriz\n" << mat << endl;
    
    Matriz<float> inv(mat);
    
    if(inv == mat) cout << "iguales\n" << endl;
    
    inv = mat.inversa();
    
    cout << "Inversa\n" << inv << endl;
    
    Matriz<float> ident(mat);
    
    ident = mat * inv;
    
    cout << "Matriz * Inversa = Identidad\n" << ident << endl;
    
    VectorMath<float> vec(2);
    vec[0] = 1;
    vec[1] = 2;
    
    cout << "Vector\n" << vec[0] << " " << vec[1] << "\n\n";
    
    VectorMath<float> vec2(2);
    
    vec2 = mat * vec;
    
    cout << "Matriz * Vector\n" << vec2[0] << " " << vec2[1] << "\n\n";
    
    return true;
}
